"""
    QualityControl

Comprehensive multi-level quality control system for atmospheric observations used in GSI.
This module implements the sophisticated quality control algorithms that ensure only 
high-quality observations are used in the data assimilation process.

# Quality Control Framework

The quality control system operates in multiple sequential stages:

1. **Gross Error Checks**: Physical bounds and consistency validation
2. **Spatial Consistency (Buddy Checks)**: Neighboring observation validation
3. **Background Checks**: Innovation-based quality assessment
4. **Variational Quality Control**: Cost function-based outlier detection
5. **Instrument-Specific Checks**: Specialized validation for different platforms

# Mathematical Foundation

## Background Check
The background check evaluates the innovation (observation minus background):
```
d = y^o - H(x^b)
```

The normalized innovation is computed as:
```
d_norm = d / sqrt(R + H*B*H^T)
```

Where:
- `y^o`: Observation value
- `H(x^b)`: Background value in observation space  
- `R`: Observation error variance
- `B`: Background error covariance
- `H`: Observation operator

## Variational Quality Control  
Uses the Anderson-Järvinen cost function modification:
```
J_QC = Σ ρ(d_i^2 / (R_i + H_i*B*H_i^T))
```

Where ρ is a robust cost function that reduces the influence of outliers.

# Quality Control Stages

Each observation progresses through multiple QC stages with increasing sophistication:
- **Stage 0**: Format validation and gross error checks
- **Stage 1**: Spatial consistency and buddy checks  
- **Stage 2**: Background departure analysis
- **Stage 3**: Variational quality control
- **Stage 4**: Final quality assessment and flagging
"""
module QualityControl

using LinearAlgebra
using Statistics  
using Distances
using ..ObservationTypes
using ..GSICoreAnalysis: AbstractAnalysisConfig, AbstractStateVector

# Export main functions and types
export QualityControlConfig, QualityControlResult, QualityControlStats
export apply_quality_control, gross_error_check, buddy_check, background_check, variational_qc
export compute_innovation_statistics, update_observation_errors, compute_qc_statistics
export BuddyCheckConfig, BackgroundCheckConfig, VariationalQCConfig

"""
    QualityControlConfig{T<:AbstractFloat}

Configuration parameters for the quality control system.

# Fields
- `gross_error_bounds::Dict{ObservationVariable,Tuple{T,T}}`: Physical bounds for each variable
- `buddy_check_radius::T`: Spatial radius for buddy checks (meters)
- `buddy_check_threshold::T`: Threshold for buddy check rejection
- `background_check_threshold::T`: Threshold for background check rejection  
- `variational_qc_threshold::T`: Threshold for variational QC rejection
- `min_observations_buddy::Int`: Minimum observations required for buddy check
- `max_background_departure::T`: Maximum allowed background departure (normalized)
- `enable_adaptive_qc::Bool`: Enable adaptive quality control thresholds
- `qc_stages::Vector{Symbol}`: Quality control stages to apply
- `platform_specific_config::Dict{String,Dict{Symbol,T}}`: Platform-specific QC parameters
"""
struct QualityControlConfig{T<:AbstractFloat}
    gross_error_bounds::Dict{ObservationVariable,Tuple{T,T}}
    buddy_check_radius::T
    buddy_check_threshold::T  
    background_check_threshold::T
    variational_qc_threshold::T
    min_observations_buddy::Int
    max_background_departure::T
    enable_adaptive_qc::Bool
    qc_stages::Vector{Symbol}
    platform_specific_config::Dict{String,Dict{Symbol,T}}
    
    function QualityControlConfig{T}() where T
        # Default physical bounds for atmospheric variables
        gross_bounds = Dict{ObservationVariable,Tuple{T,T}}(
            TEMPERATURE => (T(150.0), T(350.0)),           # 150K to 350K
            PRESSURE => (T(1000.0), T(110000.0)),          # 10 hPa to 1100 hPa
            SURFACE_PRESSURE => (T(50000.0), T(110000.0)), # 500 hPa to 1100 hPa
            WIND_U => (T(-150.0), T(150.0)),               # ±150 m/s
            WIND_V => (T(-150.0), T(150.0)),               # ±150 m/s
            SPECIFIC_HUMIDITY => (T(0.0), T(0.05)),        # 0 to 50 g/kg
            BRIGHTNESS_TEMPERATURE => (T(100.0), T(400.0)) # 100K to 400K
        )
        
        # Default platform-specific configurations
        platform_config = Dict{String,Dict{Symbol,T}}(
            "radiosonde" => Dict(:background_threshold => T(3.0), :buddy_threshold => T(2.5)),
            "aircraft" => Dict(:background_threshold => T(3.5), :buddy_threshold => T(2.0)),
            "surface" => Dict(:background_threshold => T(2.5), :buddy_threshold => T(2.0)),
            "satellite" => Dict(:background_threshold => T(4.0), :buddy_threshold => T(3.0))
        )
        
        new{T}(
            gross_bounds,
            T(100_000),    # buddy_check_radius (100 km)
            T(2.0),        # buddy_check_threshold
            T(3.0),        # background_check_threshold  
            T(4.0),        # variational_qc_threshold
            3,             # min_observations_buddy
            T(5.0),        # max_background_departure
            true,          # enable_adaptive_qc
            [:gross_error, :buddy_check, :background_check, :variational_qc],
            platform_config
        )
    end
end

# Convenience constructor
QualityControlConfig() = QualityControlConfig{Float64}()

"""
    BuddyCheckConfig{T<:AbstractFloat}

Configuration for spatial consistency (buddy) checks.

# Fields
- `search_radius::T`: Spatial search radius for finding buddy observations
- `height_tolerance::T`: Maximum height difference for vertical buddy checks
- `time_tolerance::T`: Maximum time difference for temporal buddy checks  
- `min_buddies::Int`: Minimum number of buddy observations required
- `rejection_threshold::T`: Threshold for buddy check rejection
- `distance_weighting::Bool`: Use distance-based weighting in buddy check
- `elevation_correction::Bool`: Apply elevation correction for surface observations
"""
struct BuddyCheckConfig{T<:AbstractFloat}
    search_radius::T
    height_tolerance::T
    time_tolerance::T
    min_buddies::Int
    rejection_threshold::T
    distance_weighting::Bool
    elevation_correction::Bool
    
    function BuddyCheckConfig{T}() where T
        new{T}(
            T(100_000),    # search_radius (100 km)
            T(500.0),      # height_tolerance (500 m)
            T(3600.0),     # time_tolerance (1 hour)
            3,             # min_buddies
            T(2.0),        # rejection_threshold
            true,          # distance_weighting
            true           # elevation_correction
        )
    end
end

BuddyCheckConfig() = BuddyCheckConfig{Float64}()

"""
    BackgroundCheckConfig{T<:AbstractFloat}

Configuration for background departure checks.

# Fields
- `rejection_threshold::T`: Threshold for background check rejection (normalized units)
- `warning_threshold::T`: Threshold for background check warning
- `adaptive_threshold::Bool`: Use adaptive thresholds based on local conditions
- `error_inflation_factor::T`: Factor to inflate observation errors for failed checks
- `use_ensemble_spread::Bool`: Use ensemble spread in background check
"""
struct BackgroundCheckConfig{T<:AbstractFloat}
    rejection_threshold::T
    warning_threshold::T
    adaptive_threshold::Bool
    error_inflation_factor::T
    use_ensemble_spread::Bool
    
    function BackgroundCheckConfig{T}() where T
        new{T}(
            T(3.0),        # rejection_threshold
            T(2.0),        # warning_threshold  
            true,          # adaptive_threshold
            T(2.0),        # error_inflation_factor
            false          # use_ensemble_spread
        )
    end
end

BackgroundCheckConfig() = BackgroundCheckConfig{Float64}()

"""
    VariationalQCConfig{T<:AbstractFloat}

Configuration for variational quality control.

# Fields
- `cost_threshold::T`: Cost function threshold for observation rejection
- `influence_radius::T`: Spatial influence radius for variational QC
- `max_iterations::Int`: Maximum iterations for variational QC optimization
- `convergence_tolerance::T`: Convergence tolerance for iterations
- `robust_cost_function::Symbol`: Type of robust cost function (:huber, :bisquare, :fair)
- `outlier_probability::T`: A priori outlier probability
"""
struct VariationalQCConfig{T<:AbstractFloat}
    cost_threshold::T
    influence_radius::T
    max_iterations::Int
    convergence_tolerance::T
    robust_cost_function::Symbol
    outlier_probability::T
    
    function VariationalQCConfig{T}() where T
        new{T}(
            T(4.0),           # cost_threshold
            T(200_000),       # influence_radius (200 km)
            10,               # max_iterations
            T(0.01),         # convergence_tolerance
            :huber,          # robust_cost_function
            T(0.1)           # outlier_probability (10%)
        )
    end
end

VariationalQCConfig() = VariationalQCConfig{Float64}()

"""
    QualityControlResult{T<:AbstractFloat}

Results from quality control processing.

# Fields
- `observation::AbstractObservation{T}`: Original observation with updated QC flags
- `qc_stage_results::Dict{Symbol,Bool}`: Results from each QC stage
- `innovation::T`: Innovation (observation minus background)
- `normalized_innovation::T`: Normalized innovation
- `buddy_statistics::Union{T,Nothing}`: Buddy check statistics
- `background_statistics::Union{T,Nothing}`: Background check statistics
- `variational_cost::Union{T,Nothing}`: Variational QC cost
- `final_decision::Symbol`: Final QC decision (:accept, :reject, :monitor)
- `rejection_reason::Union{String,Nothing}`: Reason for rejection if applicable
"""
struct QualityControlResult{T<:AbstractFloat}
    observation::AbstractObservation{T}
    qc_stage_results::Dict{Symbol,Bool}
    innovation::T
    normalized_innovation::T
    buddy_statistics::Union{T,Nothing}
    background_statistics::Union{T,Nothing}
    variational_cost::Union{T,Nothing}
    final_decision::Symbol
    rejection_reason::Union{String,Nothing}
end

"""
    QualityControlStats{T<:AbstractFloat}

Overall quality control statistics for a batch of observations.

# Fields
- `total_observations::Int`: Total number of observations processed
- `accepted_observations::Int`: Number of observations accepted for analysis
- `rejected_observations::Int`: Number of observations rejected
- `monitor_observations::Int`: Number of observations set to monitoring only
- `rejection_reasons::Dict{String,Int}`: Count of observations rejected by reason
- `mean_innovation::T`: Mean innovation across all observations
- `rms_innovation::T`: RMS innovation across all observations
- `bias_statistics::Dict{ObservationVariable,T}`: Bias statistics by variable type
"""
struct QualityControlStats{T<:AbstractFloat}
    total_observations::Int
    accepted_observations::Int
    rejected_observations::Int
    monitor_observations::Int
    rejection_reasons::Dict{String,Int}
    mean_innovation::T
    rms_innovation::T
    bias_statistics::Dict{ObservationVariable,T}
end

"""
    apply_quality_control(observations::Vector{AbstractObservation{T}},
                         background::AbstractStateVector{T},
                         config::QualityControlConfig{T}) where T

Apply comprehensive quality control to a collection of observations.

# Arguments
- `observations::Vector{AbstractObservation{T}}`: Input observations
- `background::AbstractStateVector{T}`: Background state for innovation computation
- `config::QualityControlConfig{T}`: Quality control configuration

# Returns
- `Vector{QualityControlResult{T}}`: QC results for each observation
- `QualityControlStats{T}`: Overall QC statistics

# Examples
```julia
config = QualityControlConfig{Float64}()
results, stats = apply_quality_control(observations, background_state, config)

# Filter for accepted observations
accepted_obs = [r.observation for r in results if r.final_decision == :accept]
```
"""
function apply_quality_control(observations::Vector{AbstractObservation{T}},
                              background::AbstractStateVector{T},
                              config::QualityControlConfig{T}) where T
    results = QualityControlResult{T}[]
    
    # Initialize tracking statistics
    total_count = length(observations)
    accepted_count = 0
    rejected_count = 0
    monitor_count = 0
    rejection_reasons = Dict{String,Int}()
    innovations = T[]
    
    # Process each observation through QC stages
    for obs in observations
        qc_result = process_single_observation(obs, observations, background, config)
        push!(results, qc_result)
        
        # Update statistics
        if qc_result.final_decision == :accept
            accepted_count += 1
        elseif qc_result.final_decision == :reject
            rejected_count += 1
            if qc_result.rejection_reason !== nothing
                rejection_reasons[qc_result.rejection_reason] = 
                    get(rejection_reasons, qc_result.rejection_reason, 0) + 1
            end
        else
            monitor_count += 1
        end
        
        push!(innovations, qc_result.innovation)
    end
    
    # Compute overall statistics
    mean_innovation = mean(innovations)
    rms_innovation = sqrt(mean(innovations.^2))
    
    # Compute bias statistics by variable type
    bias_stats = compute_bias_statistics(results)
    
    stats = QualityControlStats{T}(
        total_count, accepted_count, rejected_count, monitor_count,
        rejection_reasons, mean_innovation, rms_innovation, bias_stats
    )
    
    return results, stats
end

"""
    process_single_observation(obs::AbstractObservation{T}, 
                              all_observations::Vector{AbstractObservation{T}},
                              background::AbstractStateVector{T},
                              config::QualityControlConfig{T}) where T

Process a single observation through all quality control stages.
"""
function process_single_observation(obs::AbstractObservation{T}, 
                                   all_observations::Vector{AbstractObservation{T}},
                                   background::AbstractStateVector{T},
                                   config::QualityControlConfig{T}) where T
    qc_stage_results = Dict{Symbol,Bool}()
    final_decision = :accept
    rejection_reason = nothing
    
    # Compute innovation (observation minus background)
    innovation = compute_innovation(obs, background)
    normalized_innovation = innovation / obs.observation_error.total_error
    
    # Stage 1: Gross error check
    if :gross_error in config.qc_stages
        gross_check_passed = gross_error_check(obs, config)
        qc_stage_results[:gross_error] = gross_check_passed
        
        if !gross_check_passed
            final_decision = :reject
            rejection_reason = "gross_error_bounds"
        end
    end
    
    # Stage 2: Buddy check (only if gross error passed)
    buddy_stats = nothing
    if final_decision == :accept && :buddy_check in config.qc_stages
        buddy_check_passed, buddy_stats = buddy_check(obs, all_observations, config)
        qc_stage_results[:buddy_check] = buddy_check_passed
        
        if !buddy_check_passed
            final_decision = :reject  
            rejection_reason = "buddy_check_failed"
        end
    end
    
    # Stage 3: Background check (only if previous stages passed)
    background_stats = nothing
    if final_decision == :accept && :background_check in config.qc_stages
        background_check_passed, background_stats = background_check(obs, background, config)
        qc_stage_results[:background_check] = background_check_passed
        
        if !background_check_passed
            # Allow monitoring for moderate background check failures
            if abs(normalized_innovation) > config.background_check_threshold * 1.5
                final_decision = :reject
                rejection_reason = "background_check_failed"
            else
                final_decision = :monitor
            end
        end
    end
    
    # Stage 4: Variational QC (only if observation not rejected)
    variational_cost = nothing
    if final_decision != :reject && :variational_qc in config.qc_stages
        variational_passed, variational_cost = variational_qc(obs, background, config)
        qc_stage_results[:variational_qc] = variational_passed
        
        if !variational_passed
            final_decision = :reject
            rejection_reason = "variational_qc_failed"
        end
    end
    
    # Update observation QC flags based on results
    updated_obs = update_observation_qc_flags(obs, qc_stage_results, final_decision)
    
    return QualityControlResult{T}(
        updated_obs, qc_stage_results, innovation, normalized_innovation,
        buddy_stats, background_stats, variational_cost, 
        final_decision, rejection_reason
    )
end

"""
    gross_error_check(obs::AbstractObservation{T}, config::QualityControlConfig{T}) where T

Perform gross error check based on physical bounds.

# Arguments
- `obs::AbstractObservation{T}`: Observation to check
- `config::QualityControlConfig{T}`: QC configuration with bounds

# Returns
- `Bool`: True if observation passes gross error check

# Examples
```julia
config = QualityControlConfig{Float64}()
passed = gross_error_check(temperature_obs, config)
```
"""
function gross_error_check(obs::AbstractObservation{T}, config::QualityControlConfig{T}) where T
    # Check if bounds are defined for this variable type
    if !haskey(config.gross_error_bounds, obs.variable)
        return true  # No bounds defined, pass by default
    end
    
    min_bound, max_bound = config.gross_error_bounds[obs.variable]
    
    # Check basic physical bounds
    if obs.observed_value < min_bound || obs.observed_value > max_bound
        return false
    end
    
    # Additional checks for specific observation types
    if obs isa RadianceObservation
        # Check for satellite radiance specific bounds
        if obs.channel <= 0 || obs.scan_angle < 0 || obs.scan_angle > 90
            return false
        end
    elseif obs isa RadiosondeObservation  
        # Check for radiosonde specific bounds
        if obs.pressure_level <= 0 || obs.time_since_launch < 0
            return false
        end
    end
    
    return true
end

"""
    buddy_check(obs::AbstractObservation{T}, 
               all_observations::Vector{AbstractObservation{T}},
               config::QualityControlConfig{T}) where T

Perform spatial consistency (buddy) check against neighboring observations.

# Arguments
- `obs::AbstractObservation{T}`: Observation to check
- `all_observations::Vector{AbstractObservation{T}}`: All available observations
- `config::QualityControlConfig{T}`: QC configuration

# Returns
- `Tuple{Bool,T}`: (passed, buddy_statistic)

The buddy check compares an observation with nearby observations of the same type.
"""
function buddy_check(obs::AbstractObservation{T}, 
                    all_observations::Vector{AbstractObservation{T}},
                    config::QualityControlConfig{T}) where T
    
    # Find buddy observations within specified radius and time window
    buddy_config = BuddyCheckConfig{T}()
    buddy_observations = find_buddy_observations(obs, all_observations, buddy_config)
    
    if length(buddy_observations) < buddy_config.min_buddies
        return true, zero(T)  # Not enough buddies, pass by default
    end
    
    # Compute buddy statistic
    buddy_values = [buddy_obs.observed_value for buddy_obs in buddy_observations]
    buddy_distances = [geographic_distance(obs.location, buddy_obs.location) 
                      for buddy_obs in buddy_observations]
    
    # Distance-weighted mean of buddy observations
    if buddy_config.distance_weighting
        weights = exp.(-buddy_distances ./ (buddy_config.search_radius / 3))
        weights ./= sum(weights)
        buddy_mean = sum(weights .* buddy_values)
        buddy_std = sqrt(sum(weights .* (buddy_values .- buddy_mean).^2))
    else
        buddy_mean = mean(buddy_values)
        buddy_std = std(buddy_values)
    end
    
    # Compute normalized departure from buddy mean
    buddy_departure = abs(obs.observed_value - buddy_mean)
    
    # Account for observation and buddy uncertainties
    buddy_uncertainty = sqrt(obs.observation_error.total_error^2 + (buddy_std^2 / length(buddy_observations)))
    normalized_departure = buddy_departure / buddy_uncertainty
    
    # Apply buddy check threshold
    passed = normalized_departure <= buddy_config.rejection_threshold
    
    return passed, normalized_departure
end

"""
    background_check(obs::AbstractObservation{T}, 
                    background::AbstractStateVector{T},
                    config::QualityControlConfig{T}) where T

Perform background departure check.

# Arguments
- `obs::AbstractObservation{T}`: Observation to check
- `background::AbstractStateVector{T}`: Background state  
- `config::QualityControlConfig{T}`: QC configuration

# Returns
- `Tuple{Bool,T}`: (passed, normalized_innovation)

The background check evaluates the innovation against expected background error.
"""
function background_check(obs::AbstractObservation{T}, 
                         background::AbstractStateVector{T},
                         config::QualityControlConfig{T}) where T
    
    # Compute innovation
    innovation = compute_innovation(obs, background)
    
    # Estimate total error (observation + background)
    # Simplified - in practice would use full B matrix
    background_error = obs.observation_error.total_error  # Simplified estimate
    total_error = sqrt(obs.observation_error.total_error^2 + background_error^2)
    
    # Compute normalized innovation
    normalized_innovation = abs(innovation) / total_error
    
    # Get threshold (platform-specific if available)
    threshold = config.background_check_threshold
    platform_id = get_platform_id(obs)
    if haskey(config.platform_specific_config, platform_id)
        platform_config = config.platform_specific_config[platform_id]
        if haskey(platform_config, :background_threshold)
            threshold = platform_config[:background_threshold]
        end
    end
    
    # Apply adaptive threshold if enabled
    if config.enable_adaptive_qc
        # Adjust threshold based on local observation density and conditions
        adaptive_factor = compute_adaptive_factor(obs, config)
        threshold *= adaptive_factor
    end
    
    passed = normalized_innovation <= threshold
    
    return passed, normalized_innovation
end

"""
    variational_qc(obs::AbstractObservation{T}, 
                  background::AbstractStateVector{T},
                  config::QualityControlConfig{T}) where T

Perform variational quality control using robust cost function.

# Arguments
- `obs::AbstractObservation{T}`: Observation to check
- `background::AbstractStateVector{T}`: Background state
- `config::QualityControlConfig{T}`: QC configuration

# Returns
- `Tuple{Bool,T}`: (passed, variational_cost)

Variational QC uses a robust cost function to identify outliers in the context
of the full observation-analysis system.
"""
function variational_qc(obs::AbstractObservation{T}, 
                       background::AbstractStateVector{T},
                       config::QualityControlConfig{T}) where T
    
    # Compute innovation
    innovation = compute_innovation(obs, background)
    
    # Compute normalized squared innovation
    total_error = obs.observation_error.total_error
    normalized_squared_innovation = (innovation / total_error)^2
    
    # Apply robust cost function
    variational_config = VariationalQCConfig{T}()
    robust_cost = apply_robust_cost_function(
        normalized_squared_innovation, 
        variational_config.robust_cost_function,
        variational_config.cost_threshold
    )
    
    # Decision based on robust cost
    passed = robust_cost <= variational_config.cost_threshold
    
    return passed, robust_cost
end

# Helper functions

"""
    find_buddy_observations(obs::AbstractObservation{T}, 
                           all_observations::Vector{AbstractObservation{T}},
                           config::BuddyCheckConfig{T}) where T

Find observations suitable for buddy checking.
"""
function find_buddy_observations(obs::AbstractObservation{T}, 
                                all_observations::Vector{AbstractObservation{T}},
                                config::BuddyCheckConfig{T}) where T
    buddies = AbstractObservation{T}[]
    
    for other_obs in all_observations
        # Skip self
        if other_obs === obs
            continue
        end
        
        # Must be same variable type
        if other_obs.variable != obs.variable
            continue
        end
        
        # Check spatial distance
        distance = geographic_distance(obs.location, other_obs.location)
        if distance > config.search_radius
            continue
        end
        
        # Check time difference
        time_diff = abs(Dates.value(other_obs.observation_time - obs.observation_time)) / 1000.0  # seconds
        if time_diff > config.time_tolerance
            continue
        end
        
        # Check vertical separation (for upper-air observations)
        if obs isa RadiosondeObservation && other_obs isa RadiosondeObservation
            height_diff = abs(other_obs.location.altitude - obs.location.altitude)
            if height_diff > config.height_tolerance
                continue
            end
        end
        
        push!(buddies, other_obs)
    end
    
    return buddies
end

"""
    geographic_distance(loc1::ObservationLocation{T}, loc2::ObservationLocation{T}) where T

Compute geographic distance between two locations using haversine formula.
"""
function geographic_distance(loc1::ObservationLocation{T}, loc2::ObservationLocation{T}) where T
    # Earth radius in meters
    R = T(6371000)
    
    # Convert to radians
    lat1_rad = deg2rad(loc1.latitude)
    lat2_rad = deg2rad(loc2.latitude) 
    lon1_rad = deg2rad(loc1.longitude)
    lon2_rad = deg2rad(loc2.longitude)
    
    # Haversine formula
    dlat = lat2_rad - lat1_rad
    dlon = lon2_rad - lon1_rad
    
    a = sin(dlat/2)^2 + cos(lat1_rad) * cos(lat2_rad) * sin(dlon/2)^2
    c = 2 * atan(sqrt(a), sqrt(1-a))
    
    return R * c
end

"""
    compute_innovation(obs::AbstractObservation{T}, background::AbstractStateVector{T}) where T

Compute innovation (observation minus background) for an observation.
"""
function compute_innovation(obs::AbstractObservation{T}, background::AbstractStateVector{T}) where T
    # Simplified innovation computation
    # In practice, would use forward operator H to compute H(x_b)
    return obs.observed_value - obs.background_value
end

"""
    apply_robust_cost_function(x::T, func_type::Symbol, threshold::T) where T

Apply robust cost function for variational quality control.
"""
function apply_robust_cost_function(x::T, func_type::Symbol, threshold::T) where T
    if func_type == :huber
        return x <= threshold ? x : 2 * sqrt(threshold * x) - threshold
    elseif func_type == :bisquare
        return x <= threshold ? x : threshold
    elseif func_type == :fair
        c = sqrt(threshold)
        return c^2 * (abs(x/c) - log(1 + abs(x/c)))
    else
        return x  # Linear (no robustification)
    end
end

"""
    compute_adaptive_factor(obs::AbstractObservation{T}, config::QualityControlConfig{T}) where T

Compute adaptive factor for quality control thresholds.
"""
function compute_adaptive_factor(obs::AbstractObservation{T}, config::QualityControlConfig{T}) where T
    # Simplified adaptive factor - in practice would consider:
    # - Local observation density
    # - Data assimilation window position  
    # - Model performance in the region
    # - Seasonal/diurnal variations
    
    return T(1.0)  # No adaptation in simplified version
end

"""
    update_observation_qc_flags(obs::AbstractObservation{T}, 
                               qc_results::Dict{Symbol,Bool},
                               final_decision::Symbol) where T

Update observation quality control flags based on QC results.
"""
function update_observation_qc_flags(obs::AbstractObservation{T}, 
                                    qc_results::Dict{Symbol,Bool},
                                    final_decision::Symbol) where T
    # Create updated QC flags
    updated_flags = obs.quality_control
    
    # Update analysis usage flag
    if final_decision == :accept
        updated_flags = QualityControlFlags(
            updated_flags.gross_error_flag,
            updated_flags.background_check_flag,
            updated_flags.buddy_check_flag,
            updated_flags.domain_check_flag,
            updated_flags.time_check_flag,
            updated_flags.duplicate_flag,
            1,  # analysis_usage_flag = use in analysis
            updated_flags.observation_error_flag,
            updated_flags.preprocessing_flag,
            updated_flags.instrument_flag
        )
    elseif final_decision == :reject
        updated_flags = QualityControlFlags(
            updated_flags.gross_error_flag,
            updated_flags.background_check_flag,
            updated_flags.buddy_check_flag,
            updated_flags.domain_check_flag,
            updated_flags.time_check_flag,
            updated_flags.duplicate_flag,
            -1, # analysis_usage_flag = reject
            updated_flags.observation_error_flag,
            updated_flags.preprocessing_flag,
            updated_flags.instrument_flag
        )
    else  # monitor
        updated_flags = QualityControlFlags(
            updated_flags.gross_error_flag,
            updated_flags.background_check_flag,
            updated_flags.buddy_check_flag,
            updated_flags.domain_check_flag,
            updated_flags.time_check_flag,
            updated_flags.duplicate_flag,
            0,  # analysis_usage_flag = monitor only
            updated_flags.observation_error_flag,
            updated_flags.preprocessing_flag,
            updated_flags.instrument_flag
        )
    end
    
    # Return updated observation (would need proper constructor for each type)
    # Simplified: return original observation (in practice, would create new instance)
    return obs
end

"""
    compute_bias_statistics(results::Vector{QualityControlResult{T}}) where T

Compute bias statistics by observation variable type.
"""
function compute_bias_statistics(results::Vector{QualityControlResult{T}}) where T
    bias_stats = Dict{ObservationVariable,T}()
    
    # Group innovations by variable type
    variable_innovations = Dict{ObservationVariable,Vector{T}}()
    
    for result in results
        if result.final_decision == :accept || result.final_decision == :monitor
            var_type = result.observation.variable
            if !haskey(variable_innovations, var_type)
                variable_innovations[var_type] = T[]
            end
            push!(variable_innovations[var_type], result.innovation)
        end
    end
    
    # Compute mean bias for each variable
    for (var_type, innovations) in variable_innovations
        if length(innovations) > 0
            bias_stats[var_type] = mean(innovations)
        end
    end
    
    return bias_stats
end

"""
    get_platform_id(obs::AbstractObservation) -> String

Extract platform identifier from observation (defined in ObservationTypes).
"""
function get_platform_id(obs::SurfaceObservation)::String
    return "surface"
end

function get_platform_id(obs::RadiosondeObservation)::String
    return "radiosonde"
end

function get_platform_id(obs::AircraftObservation)::String
    return "aircraft"
end

function get_platform_id(obs::RadianceObservation)::String
    return "satellite"
end

function get_platform_id(obs::AbstractObservation)::String
    return "unknown"
end

end # module QualityControl